function DFT = fourierTransform(signal, N, Fs, out, caption)

  if (nargin < 4)
    out = false;
    caption = '';
  end
  
  % Test-case
  % t = 0:100; Fs = 1; N = length(t); f = .25;
  % signal = sin(2*pi*t*f);
  % caption = 'sin(2*pi*t*f)';
  
  % Center signal
  sig = signal-mean(signal);
  
  Ts = 1/Fs;
  Tmax = (N-1) * Ts;
  
  NFFT = N;   % 2^nextpow2(N);   % Use next power of 2 for efficiency
  DFT = fft(sig,NFFT)/NFFT;
  f = Fs/2*linspace(0,1,NFFT/2+1);
  
  % Quick Check:
  % Sum of the power spectrum should equal variance of the signal
  powsum = sum(abs(DFT).^2);
  variance = sum(sig.^2)/length(sig);
  assert(abs(powsum-variance) < 0.000000001);
  
  % Plot
  if (out)
    if ~strcmp(caption, '')
      figure('Name', caption);
    end
    subplot(2,1,1);
    plot(0:Ts:Tmax, sig);
    xlabel('Time (sec)');
    ylabel('Y(t)');
    subplot(2,1,2);
    plot(f,2*abs(DFT(1:NFFT/2+1)).^2);
    xlabel('Frequency (Hz)');
    ylabel('|Y(f)|^2');
  end

end

